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I.  INTRODUCTION  AND  BACKGROUND 


The  need  for  a definitive  formulation  of  the  relationship  between 

neutron  energy  and  displacement  damage  has  been  noted  by  Conrad.^  A 
rigorously  determined  curve  based  on  the  physics  of  the  neutron-material 
interactions  and  the  appropriate  neutron  cross  sections,  which  can  be  ex- 
tended to  other  materials  where  adequate  cross  section  information  is 
available,  is  reauired.  Experimental  evaluations  of  the  accuracy  and 
the  range  of  applicability  are  also  required.  Some  preliminary  work  on 

2 3 

this  problem  has  been  reported,  ’ and,  in  this  memorandum,  we  present  a 
more  detailed  calculation  of  the  energy  dependence  of  neutron- induced 
displacement  damage  in  silicon.  A computer  code  to  calculate  permanent 
damage,  specifically,  the  energy  available  for  displacement  damage  using 
Lindhard  theory,  has  been  prepared  and  used  to  calculate  the  damage  in 
silicon  for  incident  neutrons  from  0.1  MeV  to  20  MeV  using  a neutron 
energy  spread  of  5 keV.  Considerable  care  has  been  taken  (see  Reference 
3)  to  utilize  all  available  neutron  cross  section  information  and  the  re- 
sults presented  here  aie  based  on  silicon  cross  sections  from  a DNA  file 
tape,  DNA-4151  MOD.  3 (obtained  from  RSIC,  ORNL) . 

Early  efforts  to  calculate  the  permanent  damage  in  silicon  include 

4 

the  work  of  Smith  et  al.,  who  carefully  considered  the  nuclear  data  avail 

able  at  that  time,  and  Stein, ^ who  noted  a factor  of  two  difference  be- 
tween the  damage  calculated  for  a realistic  angular  distribution  (data  for 

aluminum  were  used)  and  that  for  an  isotropic  distribution.  Holmes^  has 


^E.E.  Conrad,  "Cnnaideration  in  Entabliahing  a Standard  for  Neutron  Dia- 
plaaement  Energy  Effeota  in  Semioonduatora , " IEEE  Trana,  Nual,  Sai.,  Vol. 
NS-18,  No.  6.  pp  200-205,  1972. 

2 

W.R.  Van  Antverp  and  J.E.  Youngblood,  "Diaplaoement  Damage  in  Silicon 
Irradiated  with  6-  to  20-MeV  Neutrona,"  Bull.  Am.  Phya.  Soa.,  Vol  27, 
p A78,  2972. 

J 

C.E.  Youngblood,  V.R.  Van  Antwerp,  and  R.M.  Tapphom,  "Diaplaoement  Dam- 
age in  Silicon  Irradtated  with  6-  to  20-MeV  Neutrona,  " BRL  Memorandum  Re- 
port MR  277,8,  April  1977. 

4 ■ 

E.C.  Smtth,  et  al. , "Theoretical  and  Experimental  Determinations  of  Neu- 
tron Energy  Deposition  in  Silicon,"  IEEE  Trans.  Nucl.  Sci,,  Vol.  NS-23, 
No.  6,  pp  12-17,  1966.  D.  Binder  et  al. , "Analytical  and  Expcrirriental 
Predictions  of  Fusion  Neutron  Radiation  Effects , " AFWL-TR-66-41,  Vol.  1, 
Hughes-Fullerton,  1966. 
b 

H.tJ.  Stein,  "Energy  Dependence  of  Neutron  Damage  in  Silicon.  " . Appl, 

Phys . , Vol.  38,  No.  2,  pp  204-210,  1967. 

R.R.  Holmes,  "Energy  Dependence  for  Carrier  Removal  and  Lifetime  Damage 
by  Fast  Neutrons  in  Silicon,"  Bell  Telephone  Laboratories  Weapons  Effects 
Studies,  Report  to  ABMDA,  Vol.  II,  Supplement  III,  pp  67-68,  October  1, 
1970. 
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made  detailed  calculations  of  damage  using  a cylindrical  shape  of  clus- 
ters of  damage  and  his  data  have  been  considered  in  the  procedures  re- 
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commended  by  Thatcher  and  Green.  ’ More  recently  Guenzer  and  Manning 

have  calculated  the  displacement  damage  in  InSb  for  1 MeV  and  14  MeV 
neutrons,  and  they  compare  these  results  with  similar  calculations  for 
silicon.  InSb  shows  very  little  energy  dependence  and  the  ratio  of  dam- 
age (14  MeV  damage/1  MeV  damage)  was  found  to  be  ca.  l.OS.  Rogers  et  al.^*^ 
have  calculated  both  ionization  and  displacement  damage  in  silicon  for 
neutron  energies  from  thermal  to  20  MeV,  These  authors  partition  the 
kinetic  energy  of  the  reaction  products  and  use  a functional  representa- 
tion of  the  Lindhard  displacement  fraction  in  calculations  similar  to 
those  presented  in  this  report.  However,  they  use  average  values  of 
emitted-particle  energy  and  the  silicon  recoil  energy  (actually  a code- 
generated average  cosine) . The  results  are  compared  with  those  of  Holmes 
and  found  to  be  about  25%  lower  in  general.  Rogers'  value  for  the  (14/1) 
ratio  is  1.8  as  compared  to  Holmes'  value  of  2.0. 

The  effectiveness  of  radiation- induced  defects  in  changing  the  elec- 
trical properties  of  a diode,  for  example,  is  certain  to  depend  on  the 
extent  of  recombination,  other  forms  of  annealing,  the  extent  and  geometry 
of  cluster  formation,  the  stability  of  the  clusters  formed,  and  the  nature 
(in  clusters  or  not)  of  the  point  defects  (divacancies,  A-centers,  etc.). 

Following  the  work  of  Gossick^^  it  is  expected  that  a space-charge  reg- 
ion and  a potential  barrier  will  form  at  a cluster  of  defects.  The  im- 
plications of  clustering  on  the  energy  dependence  of  neutron  damage  in 


7 

R.  K.  Thatcher  and  M.  L.  Green,  "TREE  Preferred  Proaeduree , " DNA  2028H, 
Edition  No,  2,  Battelle  Columbus  Laboratories,  June  2972. 

Q 

M.  L.  Green  and  R. K.  Thatcher,  "Preparation  of  a Standard  Technique  for 
Detemination  of  Neutron  Equivalence  for  Bulk  Damage  in  Silicon,  " IEEE 
Trane.  Nucl.  Sai,,  Vol.  NS-19,  No.  6,  pp  200-208,  December  1972. 

9 

C.S.  Guemer  and  Irvin  Manning,  "Calculated  Displacement  Damage  by  Neu- 
trons in  InSb,"  IEEE  Trans.  Nucl.  Sd.,  Vol.  NS-21,  No.  6,  pp  26-29, 
December  1974. 

^^V.C.  Rogers,  L,  Harris,  Jr.,  D.K.  Steinman,  and  D.E.  Bryan,  "Silicon 
Ionization  and  Displacement  Kerma  for  Neutrons  f^om  Thermal  to  20  MeV," 
IEEE  Trans.  Nucl.  Sci.,  Vol.  NS-22,  No.  6,  pp  2326-2329,  December  1975. 
(Also,  see  erratum.  Op.  Cit. , Vol.  NS-23,  No.  1,  pp  875-876,  February 
1976. 

^^B.R.  Gosaick,  "Disordered  Regions  in  Semiconductors  Bombarded  by  Fast 
Neutrons,"  J.  Appl.  Phys.,  Vol.  30,  No.  8,  pp  1214-1218,  August  1959. 
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silicon  have  recently  been  reviewed  by  van  Lint  and  Leadon.  These 

authors  comiiare  earlier  work  by  the  same  group^'^  with  available  experi- 
mental results  for  electron-,  proton-,  and  neutron- induced  radiation 
damage,  and  they  conclude  that,  when  cluster  damage  is  predominate, 
carrier  removal  is  a sublinear  function  of  recoil  energy  while  lifetime 
degradation  is  superlinear.  Specifically,  it  is  concluded  that  carrier 
removal  scales  more  like  the  total  energy  available  for  displacements 
Cor,  proportional  to  the  number  of  point  defects). 

The  evaluation  of  cluster  formation  and  annealing  has  received  con- 
tinuing interest  because  of  applications  to  ion  implantation  and  heavy- 
ion  simulation  of  radiation  damage.  Current  work  includes  that  by 

14  15 

Martv'nenko  and  by  Vinetskii  and  Kondrachuk.  The  implications  of 
both  these  works  are  that  there  is  much  diffusion  out  of  the  clusters 
that  are  formed,  that  this  results  in  a limited  significance  for  clus- 
tering, and  that  it  is  sensible  to  define  a threshold  energy  for  the 
formation  of  stable  clusters.  This  last  implication  results  from  the 
fact  that  stability  increases  with  incident  particle  energy  (defect  den- 
sity). For  neutron- induced  recoils  in  silicon  it  is  concluded^^  that 
the  stable-defect  threshold  energy  is  a silicon  recoil  of  about  5 keV. 

A complete  evaluation  of  radiation  effects,  based  on  a dynamic  model, 
can  be  done  by  computer  simulation,  as  first  demonstrated  by  Gibson  et 

al.,^^  and  an  up-to-date  review  of  these  techniques  has  been  provided 

by  Jackson  and  Morgan. This  technique  can  be  expensive  in  an  applica- 
tion to  a simple  cubic  structure  to  determine  a single  physical  parameter 

I2  ! 

V.A.J.  van  Lint  arid  R.E.  Leadon,  "In^lioations  of  Cluster  Model  of  Neu- 
tron Effects  in  Silicon,"  Lattice  Defects  in  Semiconductors,  1974," 

Inst.  Phus,  Conf.  Ser. . No.  23,  The  Institute  of  Physics,  London  and 
Brtstol,  pp  227-232,  1975. 

IZ 

V.A.J.  van  Lint,  R.E.  Leadon,  and  J.F.  Colwell,  "Energy  Dependence  of 
Displacement  Effects  in  Semiconductors,"  IEEE  Trans.  Nucl.  Sad.,  Vol. 

NS- 19,  No.  6,  pp  181-185,  1972. 

14 

Yu.  V.  Martynenko,  "Anr.e'^ling  and  Clustering  of  Defects  in  Cascades,  " 
Rad.  Effects,  Vol.  29,  pp  129-135,  1976. 

^'^V.L.  Vinetskii  and  A.V.  Kondrachuk,  "The  Kinetics  of  Formation  and  the 
Parameters  of  Radiation  Defect  Clusters  in  Silicon,"  Rad.  Effects,  Vol. 
30,  pp  227-232,  1976. 

2 S 

J.B.  Gibson,  A.N.  Goland^  M.  Milgram,  and  G.H.  Vineyard,  "Dynar.ics  of 
Radiation  Damage,"  Phys.  Rev.,  Vol.  120,  pp  1229-1253,  1960. 

17 

D.P.  Jackson  and  D.V.  Morgan  "Computer  Modelling  of  Collision  Processes 
in  Solids,"  Contemf).  Phys.,  Vol.  14,  Nc.  1,  pp  25-48,  1974. 
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several  energies,  and  the  need  to  combine  these  statistically  (the  com- 
puter simulation  itself  is  not  statistical)  limit  applications.  However, 
the  contribution  cited^  the  threshold  energy  of  displacements,  would  be 
invaluable.  An  evaluation  by  computer  simulation  of  the  'stable-defect 
threshold  energy'  suggested  above  would  also  seem  worthwhile.  Although 
potentially  useful  refinements,  these  computer  simulation  studies  are 
peripheral  to  the  present  work  and  they  have  not  been  pursued. 

A general  statement  of  the  need  for  a damage  curve  was  given  at  the 
start  of  this  report.  Some  specific  applications  are:  (a)  to  compare 
with  experiments  (as  to  confirm  results);  (b)  to  evaluate  experiments 
(for  example,  to  determine  the  energy  range  and  the  effective  average 
energy  for  damage  experiments  using  an  accelerator);  (c)  to  normalize 
simulator  and  threat  spectra  for  evaluating  simulator  experiments;  (d) 
to  evaluate  threat  spectra  in  scenarios  (for  vulnerability  analysis,  vul- 
nerability reduction,  etc.);  (e)  to  provide  a base  for  evaluation  of  the 
sensitivity  of  the  damage  curve  itself  to  the  neutron  cross  section  in- 
formation used;  and  (f)  to  evaluate  more  approximate  damage  curves  such 

1 

as  those  due  to  Messenger.^  (For  a recent  review  of  approximate  func- 
tions for  calculating  neutron-energy  equivalence,  see  Reference  8).  The 
calculations  presented  here  were  in  part  precipitated  by  a need  to  evalu- 
ate the  requirements  for  cross  section  measurements  (item  e above)  and  a 
requirement  to  correct  in  detail  for  the  spread  in  neutron  energy  encoun- 
tered in  "monoenergetic-neutron"  accelerator  radiation  effects  experiments 
(item  b above).  Even  if  interest  in  cross  section  measurements  biased 
the  procedures  chosen,  it  is  now  clear  that  the  very  detailed  treatment 
of  available  neutron  cross  section  information  has  produced  useful  informa- 
tion on  the  value  of  cross  section  details.  Also,  by  modeling  the  damage 
only  to  the  extent  of  determining  that  part  of  the  energy  available  for 
displacements,  a most-general  result  is  obtained. 


II.  COMPUTATIONAL  PROCEDURE 

Some  preliminary  calculations  on  the  energy  dependence  of  neutron 
induced  displacement  damage  in  silicon,  using  our  own  selection  of  cross 

2 3 

section  data  from  current  literature,  have  been  presented.  The  pre- 
sent report  covers  an  extension  of  those  calculations  that  includes  some 
minor  corrections,  modification  of  the  program  to  permit  direct  use  of 
any  evaluated  neutron  cross  section  tape  in  ENDF/B  format,  and  calcula- 
tions at  many  more  neutron  energies.  In  the  following  paragraphs  the 


<7.(7.  Measenger , "DiepZatyement  Damage  in  Silicon  and  Germanium  Tranais- 
tore,"  IEEE  Trma.  Nuol.  Sci.,  Vol.  NS-12,  No.  2,  pp  53-74,  April 
1965. 
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prof^ram  used  for  the  preliminary  calculations  will  be  described  and  then 
the  chanties  will  be  discussed.  The  program  presently  in  use  is  included 
as  an  apiiendix. 

A computer  program  was  prepared  that  would  accept  coefficients  for 
a Legendre  polynominal  fit  to  the  angular  distribution  of  a particular 
neutron  reaction,  determine  the  silicon  recoil  energy,  and  calculate  the 
Lindhard  fraction  of  energy  for  displacement  damage.  The  program  would 
iterate  through  60  equal  3-degree  intervals,  find  the  appropriate  solid 
angle  in  each  case,  form  the  product  of  differential  neutron  cross  sec- 
tion, recoil  energy,  displacement  fraction,  and  solid  angle,  and  sum  the 
60  contributions. 

The  calculations  of  partial  cross  section,  solid  angle,  and  silicon 
recoil  energy  are  straightforward  and  need  not  be  discussed.  The  silicon 
recoil  energy  will  be  deposited  in  the  silicon  crystal  in  both  ionization 
(inelastic  processes,  or  electronic  stopping)  and  energy  lost  to  displace 

19 

ments  (nuclear  stopping).  Lindhard,  Scharff,  and  Schiott  (LSS)  have 
suggested  that  the  ionization  scattering  process  can  be  represented  by  a 
Thomas-Fermi  potential  and  this  provides  a basis  for  determining  the  dis- 
placement fraction  of  the  recoil  energy.  To  obtain  the  displacement  frac 

20 

tion  an  approximation  was  taken  from  Bertin  et  al.  This  approximation, 
an  empirical  fit  to  the  LSS  nuclear  energy  loss,  is  given  by 


de 

dp 


1/2 

I = ^ = 

nucl  0,67  2.07e  + 0.03e' 


(1) 


where  t and  d are  the  dimensionless  parameters  defined  by  l3S  as. 


a A., 


£ = 


4tt  a N A X. 

E , p X 


Z^Z2  e^(A^  + A2)  ’ (A^  + A^)^ 


(2) 


in  terms  of  particle  energy,  E and  displacement,  X,  Here, 


19 

J.  Lindhard,  M.  Saharff,  and  H.E.  Sahiott,  "Range  Concepts  and  Heavy 
Ion  Ranges,"  Kgl.  Daneke  Videnskeb  Selskab,  Mat.  Fys.  Medd.  , Vol.  33, 
No.  14,  pp  1-42,  1963. 

20 

M.C.  Bertin,  N.  Benazer-Koller , G.G.  Seaman,  and  J.R.  MacDonald,  "Elec- 

CQ 

tronagentic  Transisition  Rates  in  Ni, " Phgs.  Rev.,  Vol.  183,  pp  964- 
977,  1969. 
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Zj,  Aj  = charj;e  and  mass  of  projectile 

Z,,  A,  = charge  and  mass  of  stopping  medium 

N = number  of  stopping  atoms  per  unit  volume 

a = 0.885.-? 

a^  = the  Bohr  radius,  0.529  x 10  ® cm 

e = the  electronic  charge 

The  LSS  electronic  energy  loss  is  taken  as  proportional  to  the  square  root 
of  the  energy  parameter, 


dc 

dp 


elec 


.1/2 


(3) 


where  the  constant  of  proportionality  as  given  by  LSS  is, 
zj''^(0.0793)(Z^  (A^  + 

The  total  energy  loss  is  the  sum  of  electronic  and  nuclear  fractions. 


dp 


tot 


de 

dp 


nucl 


elec 


(5) 


These  formulae  for  energy  loss  must  be  combined  and  integrated  to  find 
the  nuclear-stopping  fraction.  If  we  denote  by  e the  value  of  the  di- 

mensionless  energy  parameter  e for  a given  recoil  energy  E,  and  denote  by 
e some  smaller  (subsequent  in  the  scattering  sequence)  value  of  the  para- 
meter, the  fraction  of  energy  going  into  nuclear  scattering  is  given  by. 


e — I 
R dp  ^nucl 

de  I de 

'“R  e dp  'tot 

or,  by  substitution  of  (1),  (3),  and  (5)  in  (6) 


F 


(£r  - E) 


de 

1 t-  U.67k  + 2.07ke  ♦ 0.03ke‘ 


(6) 


(7) 
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1 f we  let  A = 1 + 0. b7  k 


Then , 


B = 2.07  k 
C = 0.03  k 


F 


(G,  ^ c) 


de 

A ♦ Be  ♦ Ce^ 


(8) 


The  integral  is  now  a standard  form  and  the  equation  has  the  solution, 

I 


where 


F 


1 


.2Ce  B - D.1 
'•2Ce  + B + D-'J 


D = (B^  - 4AC)^''^ 


(9) 


If  we  designate  the  argument  of  the  Zn  term  in  Equation  (9)  as  Arg 
R for  e = and  Arg  0 for  e = 0,  then  the  fraction  of  energy  for  the 

full  range  of  the  particle  is 


F 


in 


Arg  R 
Arg  0 


(10) 


Equation  (10)  is  in  a form,  with  the  constants  defined,  to  be  incorporated 
into  a computer  program,  and  F is  the  Lindhard  fraction  used  in  the  pre- 
liminary calculations. 

The  treatment  of  nuclear  stopping  as  a continuous  process  is  only 
approximate,  as  the  energy  lost  in  a single  collision  can  be  a consider- 
able fraction  of  the  total  recoil  energy.  Also,  in  range  calculations 
the  energy  losses  in  individual  collisions  are  spread,  and,  if  they  are 
treated  as  randomly  distributed  (not  the  true  distribution) , then  a 
Gaussian  distribution  of  ranges  is  obtained  for  a single  initial  recoil 
energy.  A similar  spread  must  occur  in  the  nuclear-stopping  fraction. 

At  the  same  time,  LSS  has  provided  the  fraction  of  energy  of  the  primary 
recoil  that  goes  into  nuclear  scattering  and  no  correction  was  made  in 
the  preliminary  calculations  for  multiple  collisions.  That  is,  no  cor- 
rection was  made  for  energy  lost  to  ionization  in  secondary,  tertiary, 
et  sequence,  collisions. 

This  completes  the  discussion  of  the  preliminary  calculations  and 
the  following  paragraphs  will  discuss  the  calculational  procedures  used 
currently.  The  discussion  will  center  on  changes  in  the  procedure  pre- 
sented above. 
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The  calculations  of  jiartial  cross  section,  solid  angle,  and  pri- 
mary recoil  were  done  in  the  same  way  as  the  preliminary  calculations, 
but  with  some  added  bookkeeping  for  obtaining  the  neutron  cross  sections 
and  angular  distributions  from  a DNA  tape.  The  calculation  of  the  frac- 
tion of  the  primary  recoil  energy  available  for  displacements  was  changed 
substantially  for  the  present  calculations.  A correction  was  made  for 
the  energy  lost  to  ionization  in  the  collision  sequence  and,  for  this,  a 
numerical-integration  technique  was  used  in  an  iterative  process.  As 
many  as  30  collisions  were  followed  for  the  higher-energy  recoils.  The 
technique  was  done  for  silicon  recoils  (in  silicon)  for  use  with  (n,n) 
and  (n,  n')  reactions.  Aluminum-,  magnesium-,  proton-,  and  alpha-recoils 
were  also  done  to  provide  for  (n,a)  and  (n,p)  reactions.  A more  detailed 
report  on  calculations  is  planned  by  one  of  the  authors  (JEY) . The  sig- 
nificance of  the  multiple-collision  correction  can  be  seen  by  inspection 
of  Table  I.  The  values  of  the  Lindhard  fraction  calculated  by  Equation 
(10)  above,  the  values  obtained  with  the  technique  just  described,  and 

21 

values  tabulated  by  Mayer  et  al.  are  given  for  a wide  range  of  (sili- 
con) primary  recoil  energies,  E^,  in  Table  I.  It  is  clear  that  the  mul- 
tiple collision  correction  can  be  substantial  for  high  energy  neutrons. 


TABLE  I.  Comparison  of  Estimates  of  Displacement  Fractions 
at  Several  Recoil  Energies 


E^  (keV) 

Ej  (keV)^ 

E^  (keV)‘^ 

E^  (keV)'^ 

1 

0.91 

0.79 

0.83 

3 

2.7 

2.3 

2.4 

10 

8.8 

7.0 

7.3 

30 

25. 

19. 

19. 

100 

68. 

50. 

51. 

300 

141  . 

101 . 

100. 

1000 

240. 

167. 

170. 

a.  The  energy  going  into 
Equation  (10). 

displacements,  E^,  as 

calculated  using 

b.  Displacement  energy  as  calculated  for  this  report. 

c.  Displacement  energy  as  tabulated  by  Mayer  et  al.^^ 


Mayer,  L. 
conductors, " p 


Eriksson,  and  J.A.  Davies,  "Ion  Tmfjlantation  in  Semi- 
68,  Academic  Press,  NY,  1970. 
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I 


A comparison  of  the  present  results  with  the  preliminary  calcula- 
tions'^ and  with  the  calculations  of  others  indicates  general  agreement; 
however,  at  least  four  factors  that  affect  these  comparisons  must  be 
mentioned.  First,  a mistake  in  the  preliminary  calculations  led  to  the 
(n,  n)  and  (n,  n’)  reactions  being  treated  as  (n,  p)  reactions  with  a re- 
coiling Z = 13  ion.  As  might  be  expected,  this  produced  only  a very  minor 
change  in  the  calculated  damage  (3%  at  14  MeV,  2%  at  4 MeV,  1%  at  1 MeV) . 
Second,  differences  in  the  neutron  cross  sections  used,  and  in  the  methods 
of  utilization,  produce'  variations  that  can  be  quite  irregular.  Third, 
a quite  substantial  change  is  introduced  by  including  a correction  for 
multiple  collisions  in  the  Lindhard  fraction.  This  change  is  quite  regu- 
lar with  energy  and,  when  averaged  over  a complete  range  of  recoil  ener- 
gies produced  by  a single  neutron  energy,  much  less  of  a change  results 
than  the  changes  for  individual  recoil  energies  given  in  Table  I.  Finally, 
fourth,  the  proton-  and  alpha-recoil  fractions  (and  energies  and  cross 
sections)  were  tested  in  a few  cases.  These  few  cases  confirmed  that  the 
charged-particle  products  of  the  (n,  p)  and  (n,a)  reactions  contribute 
no  more  than  1 or  2%  to  the  total  damage  and  these  contributions  were  not 
included  in  the  final  calculations. 


III.  RESULTS 

The  general  nature  of  the  results  for  0.1  to  10  MeV  can  best  be  seen 
by  inspecting  Figure  1.  This  shows  the  total  calculated  damage  as  well 
as  the  contributions  from  the  various  neutron  reactions.  Similar  results 
for  the  energy  range  10  to  20  MeV  are  shown  on  Figure  2.  Again  contribu- 
tions from  the  different  neutron  reactions  are  shown. 

The  calculations  have  been  done  for  each  5 keV  interval  from  0.1  to 
20  MeV  and  the  results  are  stored  on  disc.  To  provide  maximum  flexibil- 
ity in  obtaining  a set  of  average  damage  factors  for  any  neutron  energy 
group  structure  of  interest,  the  complete  5 keV  interval  results  are  tab- 
ulated and  average  values  for  25  keV  intervals  and  100  keV  intervals  are 
also  presented.  Table  II  presents  the  calculated  damage  at  5 keV  intervals 
for  0.1  to  4 MeV.  Table  III  continues  5 keV  intervals  for  4 to  8 MeV. 
Tables  IV  and  V are  the  5 keV  calculations  for  8 to  12  and  for  12  to  lb 
MeV,  respectively.  Tables  VI  and  VII  provide  25  keV  interval  averages  for 
0.1  to  16  MeV.  Table  VIII  has  100  keV  average  values  for  the  energy  range 
0.1  to  20.  MeV.  By  using  a combination  of  these  tables,  an  average  equal 
to  one  based  on  the  complete  set  of  5 keV  results  can  be  obtained  rapidly. 

To  compare  these  results  with  the  recent  work  by  Rogers  et  al.,  as 

corrected, average  damage  values  were  calculated  for  appropriate  parts 
of  the  22-group  structure  (the  first  9 groups  include  energies  less  than 
0.1  MeV).  Both  the  present  work  and  Rogers'  data  were  normalized  to  a 
14  MeV  damage  of  1.000.  The  results  are  shown  in  Table  IX,  and  the  over- 
all agreement  is  very  good.  The  substantial  (25%)  and  irregular  varia- 
tions are  attributed  to  the  manner  in  which  the  Lindhard  fraction  is 
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Table  Vlll.  Welative  NeuTron  l>ainaj;e  (MeV-b)  AveraRcd  Over 
0.1  MeV  interval'^  - 0.1  to  20  MeV 
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TABLE  IX.  Damage  Results  in  22-Group  Structure 


Group 

Number 

A Ej^ 
(MeV) 

Damage 
(Present 
Results) 
(MeV .mb) 

Damage 
(Present 
Results) 
Gr  22=1.00 

Damage 
(Rogers ' 
Results) ^ 
Gr  22=1.00 

Ratio 
Rogers ' / 
Present 

10 

.11-. 55 

53 

0.283 

0.255 

0.90 

11 

.55-1.11 

84 

0.449 

0.387 

0.86 

12 

1.11-1.83 

97 

0.519 

0.405 

0.78 

13 

1.83-2.35 

125 

0.668 

0.721 

1.08 

14 

2.35-2.46 

103 

0.551 

0.707 

1.28 

15 

2.46-3.01 

127 

0.679 

0.761 

1.12 

16 

3.01-4.07 

120 

0.642 

0.576 

0.90 

17 

4.07-4.97 

151 

0.807 

0.793 

0.98 

18 

4.97-6.36 

153 

0.818 

0.855 

1.05 

19 

6.36-8.19 

158 

0.845 

0.807 

0.96 

20 

8.19-10.0 

166 

0.888 

0.868 

0.98 

21 

10.0-12.2 

176 

0.941 

0.942 

1.00 

22 

12,2-15.0 

187 

1.000 

1.000 

1.00 

a 10 

Rogers,  et  al.  - all  these  values  obtained  by  dividing  published  values 
by  damage  for  Group  22,  7.36  x 10  Si  Kerma  [rad  Si/(n/cm^)]. 
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avera>;ed  tor  Uittercnt  recoil  energies.  There  is  a continuing  interest 
in  defining;  a l-MeV  etjuivalence  for  different  neutron  energies  and  sjiec- 
tra.  Several  damaj^e  ratios,  for  14  McV/1  MeV  damage,  are  shown  in  Table 
X.  The  first  four  values  indicate  how  a 1-MeV  equivalent  definition 
would  depend  on  the  energy  width  chosen  at  1 MeV.  The  second  set  of 
four  values  indicate  some  extreme  ratios,  not  likely  to  be  physically 
realized.  The  last  two  values  are  practical  energy  widths,  very  near  1 
MeV,  that  differ  in  damage  produced  by  a factor  of  two. 


TABl.n  X. 


t (1  MeV  -J 
n 


Damage  Ratios,  14  MeV/1  MeV 


Rat io , 


Damage , 
Damage , 


14  MeV 
E 


n 


1 .000  MeV  1 5 keV  1.91 

1.000  MeV  ± '25  keV  1.55 

1 .000  MeV  i 100  keV  1.99 

1.000  MeV  t 200  keV  2.15 

0.81.5  MeV  i 2.5' keV  1.05 
0.91.5  MeV  - 2.5  keV  1 .22 
1.15?  MeV  * 2.5  keV  4.79 
1.247  MeV  ♦ 2.5  keV  4.5b 

1.1  to  1.2  MeV  3.40 
0.9  to  1.0  MeV  1.70 


IV.  DISCUSSION 

Neutron  displacement  damage  in  silicon  has  been  calculated,  as  a 
function  of  neutron  energy,  in  a way  that  utilizes  all  details  of  the 
available  neutron  cross  section  information.  The  results  are  presented 
in  tabular  form  for  0.1  to  16  MeV  in  5 keV  intervals,  for  0.1  to  16  MeV 
averaged  over  25  keV  intervals,  and  for  0.1  to  20  MeV  averaged  over  100 
keV  intervals.  A comparison  with  other  current  work  has  been  presented, 
as  has  a table  of  14  MeV/1  MeV  damage  ratios.  In  the  detailed  treatment 
of  cross  sections,  the  l.indhard  fraction,  corrected  for  multiple  colli- 
sions, has  been  calculated  for  each  different  recoil  energy  and  for  the 
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appropriate  recoil  ion,  Si,  Al,  or  Mg.  For  those  reactions  producing  a 
proton  or  alpha  particle,  the  contribution  of  the  small  particle  has  not 
been  included.  These  small  particles  produce  only  a small  fraction  of 
the  displacement  damage,  and,  since  their  ranges  often  exceed  typical  de- 
vice dimensions,  it  would  be  necessary  to  calculate  the  balance  of  energy 
flux  of  these  particles  out  of  silicon  and  flux  into  the  silicon  from  its 
environment  (probably  not  silicon). 

The  calculation  of  energy  going  into  displacements  is  considered  to 

22 

be  quite  reliable.  Approximate  calculations  suggest  that  energy  depen- 
dence (ratios)  are  not  very  sensitive  to  the  details  of  LSS  at  high  ener- 

23 

gies,  while  the  experimental  confirmation  of  LSS  by  Sattler  and  Vook  is 
particularly  good  at  low  energies  for  the  displacement  fraction.  Treat- 
ment of  the  details  of  damage  sites,  point  defects  and  clusters,  has  pur- 
posely not  been  included.  If  these  details  do  not  create  differences  in 
the  annealing  of  damage  or  in  the  observable  effects  of  damage  such  as 
carrier  removal  and  carrier  lifetime,  then  they  do  not  need  to  be  in- 
cluded. To  the  extent  they  do  create  differences,  they  confirm  the  con- 

24 

jecture  of  Lohkamp  and  McKenzie,  that  there  is  no  single  energy-dependence 
relation . 

The  results  presented  have  been  used  to  plan  and  to  evaluate  accel- 
erator experiments  where  the  results  of  neutron  exposures  with  fairly 
narrow  energy-spread  have  been  analyzed.  Utility  in  a wide-range  of  ap- 
plications is  anticipated  and  these  efforts  (e.g.,  the  analysis  of  simu- 
lator spectra)  are  in  progress. 
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APPENDIX 


1.  Introduction 


The  computer  program  used  for  damage  calculations  follows  in  this 
appendix.  The  program  is  presented  as  presently  used,  making  use  of 
some  tabulated  material  and  containing  a mixture  of  FORTRAN  and  assembly 
language  for  use  on  an  SEL-86  computer.  Although  its  use  is  limited  in 
this  form,  it  is  presented  here  to  aid  in-house  users  and  to  aid  others 
in  formulating  similar  programs.  An  unusual  amount  of  redundancy  and 
error  checking  was  included  because  the  computer  was  in  a state  of  dynamic 
evaluation  at  the  time  the  program  was  prepared.  The  main  program  is 
followed  by  a program  preparation  routine  and  a number  of  sub-routines, 
and  these  are  largely  self-explanatory. 


.^3 


J.  Listing  of  Main  l’rogr;un 

CC-^mOn  ^(^00)♦S(^00)^F(200)»D<?00) 

CCmmONi/  I SS*»N/  I SSW5,  I SSW6 

CCmmON/WE  ACT/E3«S3«OfLCT«NN,AO»A(25) «M 1 » ME » M3 ,M4 »FLA6 ♦ XM I N 
dimension  TITLt(5)  ,AL(.;S)  tAHlES)  »AMAT(25t25)  ♦N0T3dO)  ♦ 

1 INT3  ( 1 0 ) ,NHT<*  ( 1 0)  . INT*»  ( 10) 

INTEGtW  Z2.Z^ 

BEAL  MJ  ,M?,M3.M4» 

DATA  E ILE5.CW1 ,TY1/»S  Z43523 1 00 1 7540^3200/ 

CAi L SETI0(EILE5.Ty1 ) 

CALL  P0ST39 
100  CCNTINUE 

CAi L SETIO(EILtS»CRl) 

BE  AD(‘='.  1 )mT,mtnUM,m1  » M2,  m3,  M4, 72  *74.  title 
1 ECBmAT (2I5.bE  10.0*2I5tbA4) 

CAI  L SE  T 10 (E ILE5. Ty 1 ) 

IE  (MT.r,E.2.0B.MT.LE.999)  GO  TO  105 

wBI TE (5, 31 ) MT 

31  ECPMAT(»  Illegal  mt  ='.I6) 
pause  mt 
GC  TO  100 
1 03  CCNT INUE 

lie  wHlTfC^.ll )Mr.MTNUM,TITLE»Ml»M2tM3*M4»Z2»Z4 
u ECOMAT (lH\ ,»bEACT10M  CODE  ** » 15 . 1 OX , • REV ISI ON» . 1 5 tZOX tSAB/ 

1 • Ml  « • ,E 1 0.5, lOX, tM2  »• .F 10.5, lOX, IM3  ■ • »F 1 0 . 5 » 1 OX . • m4  ■•,F10.S/ 

2 2bji.'72  15,  i5X, 2BX,  »ZB  =»,I5) 

call  e bacT (Zb.Mb.Z2,M2) 

call  E NTFP(mT,N13T0T,N1bT0T,N15T0T) 

REWIND  13  « rewind  1b 

call  SwITCm(mT) 
bEaD(13.21)M3,ZA,AwR,LIS»LFS 

21  format ( I 1 0,2E  1 1 .0,41  1 1 ) 

REAU(13,21)N13,T,(3,LT,N,NR3,NE3 
0»0/ I .Eft 

IE (NR3.GT.  10)  PAUSE  NR3 
113  REAO(13.22) (NHT3(I),INT3(I),I«1,NR3) 

22  FORMAT  (1  OX  ,(Sl  U ) 

WRI  TF  12)  7A,A„R,o,LIS,LFS,NE3 
12  ECRMATC  ZA  =•  ,EH.O,  lOX,  «AWR  * • , 1 PE  1 2 . 5 , 1 OX  , • Q =',0PF10.^/ 

1 • initial  state  =•, 15, lOX, 'FINAL  STATE  -•,15/ 

2 • NUMBER  OF  entries  IN  THE  CROSS-SECTION  TABLE  «',I6) 

IF (lT.GT. 1 .OR.T.NE.O.O)  wR 1 TE ( 6 , 1 3 ) T ,L T 

IJ  ECRMATC  nuclear  TEMPERATURE  « • , 1 PE  1 2 .5 , 1 OX  , « LT  «*,I5) 

WRITE  (6,14)NR3,  (NiT3(l)  ,INT3m  ,1-1,NR3) 

IB  ECRMATC  THERE  ARE»,I3,'  INTERPOLATION  RANGES • / ( 1 OX , 21 1 0 ) ) 

NTOTAL-O 
NBLK-0 
E-0  , 

S»0. 

REWIND  37 
I13-3« (NP3*2)/3 


34 


N3—2 

12b  READ<13.23)N13»  (Ed)  .5(1)  tl-ltB) 

23  F0RMAT(I10»6F11.0) 

IF  ( d3.EQ.N13)  GO  TO  120 
CALL  P0STl3dl3.Nl3.Nl3T0T) 

GO  TO  125 
128  113«I13*1 
N3-N3^3 

IF(E(3) .GE.1.E5)  GO  TO  130 
IF (N3^2.LT.NE3)  go  TO  125 
WRITE(5.33) 

33  FORMAT (•  there  ARE  NO  ENERGIES  ABOVE  10  KEV.d 
PAUSE  EMIN 

GO  TO  100 
130  N-3 

IF(E(1).GE.1.E5)  GO  TO  135 
E(1)»E(2)  » E(2)=E(3)  » E(3)«0.0 
S(1)»S(2)  » S(2)=S(3)  » S(3)»0.0 
N«2 

IF(E(I).GE.1.E5)  GO  TO  135 

E(1)»E(2)  I E(2)=0.0  I S(1)»S(2)  » S(2)«0.0 

N»  1 

IF  (Ed)  .GE.1.E5)  GO  TO  135 
PAUSE  ERR135 
135  CONTINUE 

E(1)*E(1)/1.E6  I E(2)»E (2)/1.E6  I E (3) »E (3) /I .E6 

NSKIP«n3-N»2 

LTOTAL-NE3-NSKIP 

LBLK» (LT0TAL»99) /loo 

EL«0.0  I EH=100.  I NL»Nh»1  I AL»0.0  I AH»0.0  t 1NT4A-2 

IFdSSW5.EQ.1.0R.ISS«6.EQ.l)  GO  TO  155 

READ(U.21)N1<*.ZAA,AWR4.LVT.LTT 

IF(LTT.EQ.l)  GO  TO  140 

WRITE (5.34) 

34  FORMAT (•  ONLY  LEGENDRE  REPRESENTATION  ALLOWED*) 

PAUSE  LTT 

GO  TO  100 

140  REA0(14.21 )N14«ZA4,AWR4«M*LCTtNK*NM 
AHAT«0.0 

IF (LVT.EQ.O)  GO  TO  150 

IF(LVT.EQ.1.AND.NK.EQ. ( (NM*1)*(NM*1) ) ) GO  TO  145 
WRITE(5.35)LVT.NK.NM 

35  FORMAT!*  LVT  =*.13.*  NK  «*.I4.*  NM  «*.I3) 

PAUSE  LVT 

145  NT-NM*! 

READ (14.25) ( (AMAT(I.J) . I« 1 .nT ) . J«1 .NT ) 

25  FORMAT (10X.6F11.0) 

150  READd4«21)Nl4.ZA4.AWR4tM« J.NR4.NE4 
IF (NR4.GT.10)  PAUSE  NR4 
READ (14.22) (NBT4(I) .INT4(I) . I-1.NR4) 
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wR  I Tt  (f>,  IS)  L VT  .lTT  ,LCT  tNK»NM«NE4tNR4,  (NBT4  ( I ) . INT4  ( I ) , I«1  ,NR4) 

15  ^OPMATC  I VT  =• . 15, 1 OXt *LTT  ■ • , 1 5 » 1 0 X » • LCT  -'.IS/ 

1 • NUMriEW  OF  fcNTRlFS  IN  CONVERSION  MATRIX  ■»,15»15X» 

? 'higmest  order  polynomial  used  •'*15/ 

) • nu^irer  of  entries  in  angular  distribution  table  is*, 15/ 

•4  • NUMBER  OF  interpolation  RANGES  I S • » 15/ ( I OX  , 2 1 1 0 ) ) 

AL=0.0 
Ah»0 . 0 

1 (NR4«?)  / 3*  (NK*5)  /6 

format ( I 1 0 . 1 I X ,F 1 1 .0 ,22X, 1 1 1/ ( 10X,6F1 l.O) ) 

151  REA0(14,^4)N1A,EL,NL,(AL(I),I«1,NL) 

IF  (N1<,.FQ.  1 14)  GO  to  152 

call  ROST  14 ( I 14,N14,NI4T0T) 

GC  TO  151 

152  I 1*,=  I 14»  1 ♦ (NL*5) /6 
FL=EL/ 1 ,FS 

15j  READ( 14,24)N1S,EM,nh, (AH(I) ,I«1,NH) 

IF  (NU.EQ.  114)  GO  TO  154 
call  POSTU  ( I 14,N14,N14T0T) 

GC  TO  15J 

154  I 1 4 X I 1 4 ♦ 1 ♦ ( NH *5 ) /6 
I h = FH/ 1 .E6 
N4  = 2 
NSCh4=1 
IM4A=  INT4  ( 1 ) 

\BT4A  = NBT4 ( 1 ) 

15o  IF (N.GE. 1 00)  GO  TO  160 

IF (N3*2.GE.NE3)  GO  TO  159 
IF ( I 1 3.GT.N13TOT)  GO  TO  300 
NA  = N*  I 
N » N * 3 

157  RE  AD ( 1 3,23) N1 3, (E ( I ) ,S ( I ) , I»NA,N) 

IF  ( I 1 J.EQ.Nl 3)  GO  TO  158 
call  P0ST13(I13,N13,N13T0T) 

GC  TO  157 
15t<  I13=113M 

t (T A) =t  (NA) / 1 ,t6  ; E (NA* 1 ) =E (NA» 1 ) /I ,E6  » E ( N ) »E ( N ) / 1 . t 6 

N3=N3«3 

IF (N3,2,LE ,NE 3)  GO  TC  155 
154  N*N*NE3-N3-? 

160  continue 

F=0.0 

D*0,0 

kR I Tf (6, 16 ) 

16  FORMAT (//•  I ENERGY  SIGMA  DAMAGE  ERROR*, 

* • £FF  dam  EFF  E4  EFF  F check* , 1 IX, *EL* ,BX, *EH* , 

* • INT4  N3  NS*  / ) 

DC  200  I I»1 , 100 

I«I  I 

IF ( I ,GT.N)  go  to  201 
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e:3«e  ( I ) 

S3»S(  I) 

165  IF(E3.LE.EH)  GO  TO  175 
IF (N4.GE.NE4)  GO  TO  310 
IF ( 1 14.GT.N14TOT)  GO  TO  3i0 
EL*EH  I NL=KH  I AL-0.0 
DO  170 

170  AL(I)»AH(I) 

17<e  AH*0.0 

«EA0(14f24)Nl^«EH.NHt (AH(D 
EH=EH/1 ,E6 

IF  (NU.EQ.  I 14)  60  TO  174 
CALL  POST14( I14,N14.N14T0T) 

GO  TO  172 

174  I14*I14*1*(NH*5)/6 
N4>NA* 1 

IF (N4.LE.NBT4A)  GO  TO  165 
IF (NSCH4,LT.NW4)  NSCH4»NSCM4*1 
IM4A*INT4  (NSCH4) 

NBT4A»NBT4 (NSCH4) 

GO  TO  165 

17b  IF(EL.LE.E3)  GO  TO  180 
WRITE (5.36) 

36  F0RMATCE3  IS  BELOW  E4  LOW'/ 

* »C  TO  TERMINATE.  R TO  IMMEDIATELY  GET  NEXT  REACTION  CODE.') 
PAUSE  E3L0W 
180  CONTINUE 

NN»MAX0 (NL .NH) 

DO  185  1*1. NN 

CALL  INTERP(E3.£L.EH.AL(I) .AH(I) .AA.INT4A) 

A(I)=AA*(I.I*1)*0.5 
185  CONTINUE 
A0  = 0.5 

CALL  FSIMPO(X) 

CALL  DAMCAL ( damage. ERROR. A VEDAM.AVEE4. FRACTION. TEST) 

I«I  I 

D ( I ) -DAMAGE 
F (I ) -FRACTION 
NUMBER-NTOTAL^  I 


N3A-NUMBER*NSKIP 

CALL  SSWTCH( 10. ISSW)  I iF ( ISSW.EQ. 1 ) GO  TO  200 
WRITE (6. 17) number. E ( I ) . S ( I ) . D ( I ) .ERROM.AVEOAM.AVEE4.F (I) 
* E L . EH . I NT4 A . N3 A . N4 


TEST. 


17  FORMAT (IX. I4.2F1 1.5. IPE 14.5, 0PF8.3.1P3E 11.2. El 0.1.5X.0P2F 1 0,4. 315) 

200  CONTINUE 

201  CONTINUE 
NN-MINO (N. 100) 

IF(NN.LE.O)  go  to  230 
NA-NTOTAl*! 

NAA-NTOTAL^NN 
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NPl  K =NHl  K ♦ I 
JPL  K=NHLK 

*iR  I TF  ( ^9)  MT  tMlNUM.NBLK  »LBLK  ,NN.LT0TAL  »NAtNAAf  Ml  ♦ M2  t M3  t m<,  , Zf?  » ZA  « U * 

* TlTLt.(F(I).l  = l«100)»(S(I)»I»l,100)f(L)(l)»I*ltlOO)t 

* (F  ( I ) , I = I , 1 00)  , JiLK 

«WITE(^.?02)NA,NAA,E(l)*S(l)«D(l)tF(l)tE(NN)tS(NN)tD(NN).F(NN) 

20t  FORMAT (/// IX. 1 30 1 IH*) / ' table  OUTPUT  FOR  N«».15»»  T0*»I5/ 

•?0x, 'FIRST  DATA  E = * . F 1 0 . 5 . 5X . • S IGMA* ' »F 1 0 . 5 . 5X . • DAMAGF- • ♦ 1 PE  1 2 . 5 , 

* 5X, 'EFFECTIVE  F = ' . E 1 2 . 5/20 X ♦ • LAST  DATA  E» • . OPF 1 0 .5 »5X . 

* ' ciGmax ' ,F 10. 5.5X. 'DAMAGE-' . 1PE12.5.5X. ‘EFFECTIVE  F»'.E12.5/ 

* 1 X . 1 30 ( iH* ) ) 

MOTAl.  xNiTOTAl  *NN 
N=N-1 00 

NA  = '<*  I 

IF  (NSk  IP.MOTAl.GE  .ME3)  60  TO  250 
IF  (N)  £;30,2?0.20S 
20-  DC  210  1 = 1. h 
E ( I ) =E ( I ♦ 1 00 ) 

210  S ( I ) =S  n ♦ 100 ) 

220  continue 

DC  215  1=NA.200 
E ( 1 ) =0.0 
215  S( I) =0.0 
GO  TO  155 

230  NN=NSn1P.NTOTAl 

IF (NN.EO.NE3)  GO  TO  250 
WRITE (5.231 )NN.NE3 

231  format ( 'POSSIBLE  ERROR  IN  ThE  NUMBER  OF  ENERGIES  CALCULATED'/ 

* ' (nSn IP.nCAlC)  =',I5.'  while  NE3  ='fI5/ 

* 'C  TO  ACCEPT  AS  IS.  R TO  RETRY  ENTIRE  REACTION') 

PAUSE  NERR 

250  continue 
N=  1 

25'  continue 

L5LX=LRLk-NHL« 

L TOT A|  XL  TOT AL-NTOTAL 
ntotai xnf  3-ntct al-nSk ip 

mj=m1 .m2-m3-mr 
1Z  = /2-.’4 

wR I TF (^ . 37) MT .MTNUM.LBLK,LT0TAL.NT0TAL*M1 t IZ 

3/  FCRmat(//'  •*  TFST  **'.5I5»F10.3fI5) 

GO  TO  100 

30i’  wRI  TF  (5,301  )E  (N) 

301  FORMAT ( 'ATTEMPT  TO  READ  MORE  THAN  N13TOT.'/ 

* 'F(N)  =',1PE12.5/ 

* 'C  TO  USE  AS  IS,  R TO  RETRY  ENTIRE  REACTION.') 

PAUSF  N13 

GO  TO  160 

310  wRI TF  ( 31  1 )F  (N) 

311  FORMAT ( 'ATTF MPT  TO  READ  MORE  THAN  NIrTOT.'/ 
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* 'e(N)  ■•tlPEl2.5/ 

*'C  TO  EXTRAPOLATE.  R TO  RETRY  ENTIRE  REACTION.') 
PAUSE  NU 
GO  TO  175 
END 
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' ■ * ‘ L'li  of  -Siihroutincs 

SLiHtVOuT  iNe.  SwITCH(mT) 
common/ 1 SSWN/ I SSWS, I SSw6 

CAl  L SSWTCH (5* I 5SWS)  » CALL  SSWTCH (6 f I SSW6) 

IF ( 1SSW5.FQ. 1 .OH. ISSW6.LQ.1)  GO  TO  10 

IF(mt.fQ.16)  ISS»^6=1 

IF<mT,FQ.^?»  ISSW6=1 

1F(MT.FQ.?0)  ISS«6=1 

IF  (MT.FQ.  104)  ISS*^5=l 

IF  (MT.FO.  10?)  I5S*<6=1 

10  continue 

IF ( lSSrt5.EQ. 1 ) WRITE(6.1) 

1 FCBmaK/'  reaction  is  forced  to  be  ISOTROPIC'/) 

IF  ( ISSW6.EQ.  1 ) »(RlTE(6t2) 

FCHMAK/'  HEAVY  PARTICLE  ASSUMED  TO  REMAIN  AT  REST  IN  C.O.M.*/) 
IF ( ISSwS.NE . 1 . and. iSSwe.NE. 1)  WRITE  <6»3) 

J FCRmaTI/'  reaction  treated  NORMALLY'/) 

RE  TURN 
END 
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SUBROUTINE  SETIO  (FILE,  DEVICE) 

[This  is  an  assembly  language  routine  which  dynamically  assigns  a 
symbolic  unit  number  (FILE)  used  in  FORTRAN  input/output  statements  to  a 
physical  device  (DEVICE).  It  is  used  to  permit  use  of  descriptor  5 for 
both  card  reader  and  teletype.  This  was  a temporary  local  requirement 
because  the  number  of  symbolic  input/output  devices  was  less  than  the 
number  of  hardware  devices.] 


(SUBROUTINE  ENTER  (NITT,  N13T0T,  N14T0T,  N15T0T) 

[This  is  an  assembly  language  routine  which  reads  a BCD  tape  and 
transfers  the  card  images  for  the  desired  reaction  code  (MTT)  to  files  13 
and  14.  File  13  is  the  cross-section  data  from  DNA  file  3 and  file  14 
is  the  angular  distribution  data  from  DNA  file  4.  An  integer  card  coun- 
ter (110)  precedes  the  card  images  (20A4) . This  routine  was  a local  re- 
quirement because  the  computer  used  could  not  read  BCD  coded  magnetic 
tapes  by  means  of  its  FORTRAN  accessible  device  handlers.  N13T0T  is  the 
number  of  cards  in  file  13  (DNA  file  3);  N14T0T  is  the  number  in  file  14 
(DNA  file  4);  and  N15T0T  is  set  equal  to  zero.] 


♦ V*  , r. 
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SUHPOUTlNf  P0ST13( I WAN T t NRt AD » NTOTAL ) 

WWI TF (S.?) IwANTtNRFADfNTOTAL 

<f  FCPVArC  POST  13  WANT».Ib»»  JUST  RtAD*»IS.*  OF  ThE'.IS* 

* • RFCORDS') 

IF ( IWANT.Lt.O.OR. IwANT.GT.NTOTAL)  GO  TO  SO 
N* I wANT-NRF AD 
IF (N) 10. 15,20 
1 U N » - N 

IF (N.GT . lO.OR.  I wANT.LT.  10)  GO  TO  30 
N = N ♦ 1 

DC  12  1*1. N 
12  HACKSPACF  13 
RE  TuPN 

15  backspace  13 
2C  RETURN 
30  REWIND  13 

IF ( IWANT.EQ,  1 ) RETURN 
N* I wANT-1 
DC  35  1=1. N 
3b  READ! 13, 1 ) J 
I FCRMAT(IIO) 

IF(J.EQ.N)  RETURN 
WRITE (5.3) 

3 FORMAT!'  error  POST  13  AUTOMATIC  RETRY*) 

N* I wAnT- J 
IF (N) 10. 15,20 
5(.  CCNTiNUe 

IF ( IwanT.LE.O)  go  to  60 
WRITE (5, SI ) 

51  FCRmatc  REQUESTED  RECORD  EXCEEDS  SIZE  OF  FILE  13*/ 

* ' C TO  RETRY  R TO  RESTART') 

GC  TO  70 

6f'  WfiITF(5,61) 

61  FCRMATC  NON-ROSITIVE  RECORD  NUMBER  REQUESTED  FOR  FILE  13*/ 

* • C TO  RETRY  R TO  RESTART') 

70  PAUSE  N13  ) RETURN 

E NO 
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SURROUTINE  POS  T 1 4 ( I W AM  t NRt  AD  » NTOT  AL  ) 

WRI  TE  (S.?)  IWAMtNREADfMOTAL 

I FORMATl’  POST  14  WANT'.IS**  JUST  READ'tlS**  OF  THE*. IS* 

• • RECORDS') 

IF (iwant.le.o.or.irant.gt.ntotal)  go  to  50 

N* I WAnT-NREAD 
IF (N) 10.15.20 
10  N»-N 

IF (N.GT. lO.OR. IWANT.LT.  10)  GO  TO  30 
N»N*  1 

DO  12  I«1.N 
12  BACKSPACE  14 
RETURN 

15  BACKSPACE  14 
20  RETURN 
30  REWIND  14 

IF ( IWANT.EQ. 1 ) RETURN 
N» I WANT-1 
DC  35  I«1.N 
3b  READ(14,1)J 
1 FOPMAT(IIO) 

IF(J.EQ.N)  RETURN 
WRITE (5.3) 

3 FCRMATC  ERROR  POST  14  AUTOMATIC  RETRY*) 

N= I want- J 
IF (N) 10.15.20 

50  CONTINUE 

IF  ( IWANT.LE.O)  GO  TO  60 
WRITE (5.51 ) 

51  FORMAT!*  REQUESTED  RECORD  EXCEEDS  SIZE  OF  FILE  14*/ 

• * C TO  RETRY  R TO  RESTART*) 

GO  TO  70 

60  WRITE(5.61) 

61  FORMAT!*  NON-POSITIVE  RECORD  NUMBER  REQUESTED  FOR  FILE  14*/ 

• * C TO  RETRY  R TO  RESTART') 

70  PAUSE  N14  * RETURN 

END 


subroutine  POST39 

DATA  1FLA6/0/ 

WRITE (S.l) 

1 FORMATC  POSITION  UNIT  39  FOLLOWING  mT  . VERS  I ON  ( 2 1 3 ) • ) 

IF ( IFLAG.EQ.C)  wRITE(5f2) 

2 FORMAT!*  USE  HT"3C0  TO  REWIND  UNIT  39*/ 

* • USE  MT-0  IF  UNIT  39  IS  POSITIONED  CORRECTLY*) 

IFLAGsIO 

READ(S»3)MTA.MTNUMA 

3 F0RMAT(2I3) 

IF(MTA)5»I0.15 

b REWIND  39 
WRI TE (5.6) 

b FORMAT!*  UNIT  39  HAS  BEEN  REWOUND*) 

GO  TO  20 

10  WRITE !5. 1 1 ) 

11  FORMAT!*  UNIT  39  HAS  NOT  BEEN  MOVED*) 

GO  TO  20 

lb  REWIND  39 

lb  READ!39)mT.MTNUM.IBLK»NBLK 

IF (MT.NE.MTA.CR.MTNUM.NE.MTnUMA.OR. IRLK.NE.NBLK)  60  TO  1« 
WRI TE (5. 16) 

lb  FORMAT!*  UNIT  39  HAS  BEEN  POSITIONED  AS  REQUESTED*) 

20  CONTINUE 
wRI TE !5.21) 

21  FORMAT!*  C TO  PROCEED  R To  RE-ENTER  P0ST39*) 

PAUSE  39 

RETURN 

END 
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o o n r>  o 


SUflPOUTiNE  INTERP(E»ELtEH»SL»SH,S» INT) 

INT  = 1 POP  y CONSTANT  FROM  £L  TO  EH 

INT^a  FOR  Y = A ♦ p * X 

INT  * 3 FOR  y = A ♦ B * LOG(X) 

INT  * h FOR  LOG(Y)  » A ♦ B * X 
INT  » 5 FOR  LOG(V)  = A ♦ B * LOG(X) 

IF (EL.EQ.EH)  60  TO  10 
IF (E.EQ.EL.OR.E.EO.EH)  GO  TO  10 
IF  ( INT.LT. 1 .OR, INT.GT.5)  GO  TO  11 
5 GO  TO  (10.11«12»13,14) tlNT 

10  S»SL 

IF(E.EQ,EH)  S=SH 
RETURN 

11  S»SL* (SH-SL)*(E-EL)/(Em-EL> 

RETURN 

12  IF (E.LE.O.O.OR.EL.LE.O.O.OR.EH.LE.O.O)  go  TO  11 
S»SLtJSH-SL) *AL0G(E /EL)/ ALOG( EH/EL) 

RETURN 

13  IF (SL.LE.O.O.OR.SH.LE.0.0)  GO  TO  11 
S»AL0G(SH/SL)*(E-EL)/(EH-EL) 

S«EXP(ALOG(SL) *5) 

RETURN 

14  IF  (E.LE.O.O.OR.EL.LE.O.O.OR.EH.LE.0.0)  go  to  11 
IF (SL.LE.O.O.CR.SH.LE.0.0)  go  to  11 
S»ALOG(SH/SL) *ALOG( E/EL) /ALOG( EH/EL) 
S*EXP(AL0G(SL) *S) 

RETURN 

END 
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SLPWOUTINE  FWACT ( ZP.AP.ZTfAj) 

COMMON / T ARlF S/NT Ai.AEPS.E tab (500) «F TAB (500) 

INTEGER  ZP»ZPA,ZTtZTA 
DATA  ZPA,ZTA/1000» 1000/ 

ZZ'ELOAT (ZP) •* (2. /I. ) ♦FLOAT (ZT) ** (2,/3. ) 

AEPS«325W.»AT/  ( ( AP*  AT  ) *FLOAT  ( ZP*ZT ) *SORT  ( ZZ ) ) 
AKr0,0793*(FLOAT(ZP)»*(2./3.) ) *SQRT (FLOAT ( ZT ) ) * (AP*AT) **1 .5/ 
• (72**0. 75*AP*»1 ,5*SQRT  (AT) ) 

WRITE (6.3) ZP.ZT.AP.AT.ZZ.AEPS.AK 
3 format (2 ( I8.2A) »2F 10. 5. Fl0.5*F10.6fF10. 7) 

IF (ZP.EQ.ZPA)  RETURN 
WRITE (6,1 ) 

J FOPMATC  READ  NEW  DAMAGE  FILE') 

15  REWIND  37 

20  READ (37,END=30) ZPA.APA.ZTA, ATA.NTAB, (ETAB( I ) , I»1 ,500) ♦ 
*(FTAi(I) ,1«1,500) 

IF (ZP.NE.ZPA)  GO  TO  20 
RETURN 

30  WRITE(5,31) 

31  FORMAT!'  READ  EOF'/'  C TO  REWIND  AND  RE-READ') 

CALL  STATUS(7, ISTAT) 

PAUSE  EOF 
GO  TO  15 
END 
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SUHWOUT INE  DAMCAL (DAMAGE*ERR0R.AVE0AM»AVEE4,AVEFRACT»TEST) 
CCmmON/EXTRA/SUMI »SUM2 
DOUBLE  PRECISION  SUMl»SUM2 

CCMMON/REACT/E  tSIGMA.0.LCT»NLtP0tP(2S)  »Ml  »M?«M3.M4fELAr,«  xL  IM 


10 


20 


DAMAGE»ERROR=AVEOAm=AVEEA«AvEFRACT«TEST»0.0 
IF (SIGMA. LE. 0.0)  RETURN 
IF (flag.lt.-o.S)  go  to  10 
IF (FLAG. GT. 0.5)  GO  TO  20 

PCTMIN*0.1  1 MINPTS*100  I MAXPTSsIOOO 
XLOWER»-1.0  I XUPPER*1.0 

call  simp  (XLOi»EW*XuPPER«MINPTS.mAXPTS»PCTMIN.  ANS*ERR0R) 

AVEDAmsANS  ( AVEE<»  = SUM1  J AVEFRACTsANS/SUMl 
SUM2=SUM2-1.D0  » TEST*SUm2 
DAMAGE»SIGMA*ANS 

there  are  two  factors  of  2(PI)  WHICH  CANCEL. 

ONE  IS  THE  integral  OVER  THE  POLAR  ANGLE. 

THE  other  arises  in  THE  DEFINITION  OF  THE  LEGENDRE  EXPANSION 
USED  IN  THE  ENDF/B  FILES. 

return 

DAMAGE*©. 0 

THE  REACTION  IS  ENERGETICALLY  IMPOSSIBLE. 

erpor=ioo. 

RETURN 

DAmaGE=0.0 

THE  LIGHT  particle  IS  UOUBLE-V ALUED  IN  ENERGY. 
calculation  not  SET  UP  FOR  THIS  CONDITION 
ERR0R=I 00. 

RETURN 

END 


SOH^^OUT  INE  SIHP  (XA,X8»NL.NH,PCTMIN»ANS»EWR) 

DCUHLE  PRECISION  SUM 1 , SUM2 • SO  1 » S02 

ccmmon/extra/sumi »SUM2 

DOUBLE  PRECISION  SO»S1»X»OX 

XL  = XA 

Xh  = XR 

NMIN»NL 

NMAX=NH 

IF  (NMAX.LE  . 0)  NMAX  = 9R9<^99 
PCT=PCTMIN 

IF (PCT.LE.O.O)  PCT=1, 

ANSA=0 . 

SUM1=SUM?*0,U0 

S0=FS1MP<XL)*FSIMP(XH) 

S01=SUMI  ; S02=SUM2I  SUMI«0.D0  « SUM2*0.D0 

OX=XH-XL 
N = 2 

IF (DX.GT.0.0)  GO  TO  10 
IF (DX.EO.0.0)  GO  TO  200 

err=xl 

XL  = Xh 

xh=err 

DX  = -OX 
1C  CONTINUE 

X=XL*0.500*DX 
S1=0.D0 
20  CONTINUE 
XX  = X 

S1=S1*FSIMP(XX) 

N = N>  1 
X=X*DX 

IF (X.LT.XH)  GO  TO  ?0 
ANS= (SO *4.00*5 1 )*DX/6.D0 
IF ( ANS.EQ.O. ) GO  TO  150 
ERR=100.* ( ANSA/ANS-1 . ) 

IF (ERR.GT. 100. ) ERR=100. 

IF (ERR.LT.-IOO. ) ERR=-100. 

IF (N.L  T.NMIN)  go  to  50 
IF (N.GT.NMAX)  go  to  100 
IF (ABS(ERR) .lT.PCT)  GO  TO  lOO 
50  CONTINUE 
ANSA=ANS 
S0=S0*S1 *S1 

S01=S01*5UMl*SUMl  I S02=S02*SUM2*SUM2 
SOMl=SUM2*O.DO 
DX*0.5D0*OX 
60  TO  10 
100  CONTINUE 

SOMl* (S01*4,D0*SUM1 )*OX/6.O0l  SUM2« ( S02*4 . D0*SUM2 ) *UX/6 . 00 

if(x«.ge.xa)  return 
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% 

ANS«-ANS 

SUM1«-SUMI  I SUM?*-SUM2 
RETURN 
150  CONTINUE 

IF (N.LT.NMIN)  go  to  50 

200  continue 
ANS«0.0 
ERP«0, 

SUH1>SUM2*0.D0 

RETURN 

ENO 


function  fsimP(x) 

CCmmONi/  I SSWN/  I SSWS.  1 SSW6 

COMMON/ tables/ NT Ai,AEFS*ETAb( 500)  »FTAB(50(}) 

C0MM0N/PEACT/E»5IGmA,0.LCT  tNL»P0tP(25)  »Ml  » M2  « M3  » M4  t FL  AC,  t XL  I M 

CCMM0N/PEACT2/ET.A,BtC.U.AC,DB 

DOUBLE  PRECISION  SUMl,SUM2 

COMMON /EXTRA/ SUM  1 .SUM2 

dimension  P0LT(25) 

REAL  Ml,M2»M3tM4 
E4=m4«mi*c/(m1»m2)**2 
IF ( ISSW6.EQ. 1 ) GO  TO  17 
if(lct.ne.I)  go  TO  10 

C LCT  =1  X IS  THE  COSINE  OF  THE  LABORATORY  ANGLE 

CON=OB- ( 1 .-X*X| 

IF (CON. LT. 0.0)  GO  TO  100 
CON=SQflT (OON) 

E3  = B»ET'*(X*C0N)*»2 
GO  TO  16 
10  CONTINUE 

C LCT  = ? X IS  THE  COSINE  OF  THE  CENTER  OF  MASS  ANGLE 

E3=ET*(H*D*(AC*AC)*X) 

16  continue 

IF (E3.lt. 0.0)  GO  TO  100 
EA=ET-F3 

IF (EA.LT. 0.0)  GO  TO  100 
1/  SIG=P0 

IF  ( ISSW5.E(2.  1 .OR.  ISSX6.EQ.  1 ) GO  TO  18 

u*x 

call  LEGPOL (NL.UtPOLY) 

DO  15  1=1 »NL 
ItJ  SIG=5IG*P(I)*F0LY(I) 

1«  CONTINUE 
EPS  = E'4*AEPS 
F=FTAB( 1 ) 

IF (NTAB.lE.  1 ) GO  TO  50 
F = 1 . 

IF (EA.lt.  1 .E-20)  GO  TO  50 
V=(AlOG10(EPS)-ALOG10(ETAB(2) ) )*3l.*2. 

N*  V 

IF(N.LT.2)  N=2 
IF (N.GT.NTAB-1  ) N = NTAB-1 
20  IF (EPS.GT.ETAB (N*l ) ) GO  TO  25 
IF (EPS.lt. ETAB(N) ) GO  TO  30 

F=FTAH(N) ♦ (FTA8(N*1 )-FTAB(N) ) * (EPS-ET AB (N) ) / (ET A8 ( N* 1 ) -F T A8 (N) ) 
GO  TO  50 

25  N»N*1 

IF (N.LT.NTAB-1 ) GO  TO  20 
WRITE (6«26)EPS.£TAH(NTAB) 

26  FORMATC  REDUCED  ENERGY  EXCEEDS  TABLE.  EPS  »'»1PE12.5. 

* • while  MAXIMUM  VALUE  IN  TABLF  IS*»E12.5/) 


50 


F-FTAH (NTAi) •ETAB (nTAB> /EPS 
GO  TO  50 
30  N»N-1 

IF(N.GT.I)  go  to  20 

F»l,-(l.-FTAB(2) >*( (EPS/ETAB(2) ) •*0.166667) 
GO  TO  50 
50  CONTINUE 

SUMlsSUMl ♦E4*SI6  I SUM2*SuM2^SI6 
FSIMP«E4*F»SIG 
RETURN 
100  CONTINUE 
FSlMP-0,0 
RETURN 
END 
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n n r>  r>  o 


SUHWOUTINF.  ri>lMP0(X) 

CCMMON/Wf  ACT/E  .SIGmA.Q.LCT  «NLiP0tP  (2S)  tMl  *M2,M3tM<».FLAG*  XL  IM 
CCMMON/REACT2/ET* AtB.C.Ot AC.DB 
SEAL  m1*m2,MJ*MA 
ET=E*Q 

IF (ET.LE.0.0)  go  to  1000 
U= 1 ./ < (Ml *M2) • (m3*M4) ) 

A=m1*m4*E*u/ET 
0=m1*m3*e*U/ET 
U«M2*u*  ( 1 . ♦Ml*(5/  (m?*ET  ) ) 

IF (U.LT.0.0)  GO  TO  1000 
C=M3*U 
D=M4*U 
AC=SQRT ( A*C) 

DB=0/H 
XL IM  = -1 .0 

flag=o.o 

IF (D.LT.B.ANO.LCT.EO. 1 ) GO  TO  1010 
RETURN 

1000  continue 

C REACTION  BELO**  THRESHOLD. 

FlAG=-1  . 

return 

1010  CONTINUE 

the  enfr&t  of  The  light  particle  is  double-valued. 

THE  MINIMUM  value  OF  THE  COSINE  (XLIM)  IS  RETURNED  TO  THE  CALLING 
program.  ATTEMPTED  USE  WILL  YIELD  A CALCULATION  FOP  THE 
higher  energy  BRANCH  FOR  THE  LIGHTER  PARTICLE  (THE  LOWER 
ENERGY  branch  OF  THE  HEAVY  PARTICLE  OR  INTEREST). 

XL IM=0.0 
FLAG=1 .0 
U=1 .-UB 

IF(U.GT.O.O)  ALIM=S0RT(U) 

RETURN 
END 
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SOHROUTINE  LEGPOL (N»X.P0LY) 

DIMENSION  POLY(l) 

double  PRECISION  Y.U»V 
V«X 

IF(N)  10,20.30 
10  mRITE(6,1> 

1 FORMAT (•  ERROR  IN  CALCULATION  OF  LEGENDRE  POLYNOMIALS') 
RETURN 

20  P0LY(l)»0.  » RETURN 

30  DC  5 I-l.N 
b PCLY(I)«0. 

P0LY(1)«Y  I IFIN.EO.I)  return  I U-l.DO  I V»Y  I 
DO  40  1*1, J 

Y»ni*I*l)«X*V-I*U)/(I*l)  » U»V  » V*Y 

40  P0LY(I*1)*Y 
RETURN 
END 


I 
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